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NUMERICAL  EXPERIENCE  WITH  A  SUPERFAST 
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DeKalb,  EL  60115 


William  B.  Gragg 
Department  of  Mathematics 
Naval  Postgraduate  School 
Monterey,  CA  93943 


Abstract 


We  briefly  describe  the  Generalized  Schur  Algorithm  for  the  superfast  solution  of 
positive  definite  Toeplitz  systems  of  equations  and  its  relationship  with  Sclnir’s 
algorithm  and  the  Szego  recursions.  We  then  present  some  experimental  results 
obtained  with  our  FORTRAN  implementation  of  this  superfast  Toeplitz  solver. 
We  will  see  that  the  algorithm  displays  favorable  behavior  in  that  the  growth 
rates  of  the  resulting  residuals  and  errors  are  comparable  with  those  of  the  Szego 
recursions. 


1.  Introduction 

The  problem  of  efficiently  solving  a  linear  system  of  equations  M  x  =  b  where 
M  is  an  Hermitian  positive  definite  Toeplitz  matrix  is  of  considerable  interest  in 
many  areas  of  mathematics  and  engineering.  While  standard  Gauss  or  Choleski 
methods  for  solving  such  a  system  of  order  n  require  0(n3)  arithmetic  operations, 
there  are  several  algorithms  for  obtaining  the  solution  in  0(n'2)  operations 
(21,24,5,18,13).  We  refer  to  these  algorithms  as  fast  Toeplitz  solvers.  More 
recently,  algorithms  have  been  presented  for  solving  a  Toeplitz  system  of  equa¬ 
tions  using  0(nlog2n)  arithmetic  operations  [11,6,7,15,17,22].  These  superfast 

f  Presented  at  the  International  Conference  on  Linear  Algebra  and  Its  Applications, 
Valencia,  Spain,  September,  1987.  This  research  was  supported  in  part  by  the  National 
Science  Foundation  under  grant  DMS-8704196. 
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Toeplitz  solvers  are  all  based  on  certain  doubling  procedures.  While  they  are  more 
complicated  to  implement  than  the  fast  methods,  they  require  substantially  less 
computational  work  for  sufficiently  large  problems. 

Current  fast  and  superfast  Toeplitz  solvers  are  numerically  unstable,  and 
therefore  unreliable,  when  applied  to  arbitrary  Toeplitz  systems  [8].  However, 
Cybenko  (9,  10]  has  shown  that  the  algorithms  of  Levinson,  Durbin,  and  Trench 
are  numerically  stable  for  the  class  of  positive  definite  Toeplitz  matrices  in  the 
sense  that  accurate  results  can  be  expected  when  the  problem  is  well  conditioned. 
Stability  of  some  superfast  algorithms  can  also  be  expected  for  the  class  of  posi¬ 
tive  definite  matrices  [8].  We  remark  that  some  fast  and  superfast  algorithms  are 
also  potentially  forward  stable  for  other  classes  of  matrices,  for  example,  for 
totally  positive  Toeplitz  matrices. 

The  numerical  reliability  of  fast  and  superfast  Toeplitz  solvers  is  a  topic  of 
ongoing  research.  Since  they  are  relatively  new  on  the  scene  and,  moreover, 
implementations  of  these  algorithms  are  less  tractable  than  those  of  the  fast  algo¬ 
rithms,  numerical  results  for  superfast  Toeplitz  solvers  are  currently  rare. 
Knowledge  of  the  numerical  reliability  is  important  because  superfast  methods 
allow'  for  the  solution  of  problems  which  would  require  a  prohibitive  amount  of 
time  for  a  fast  method.  We  are  not  aware  of  any  experimental  results  for  any 
superfast  Toeplitz  solvers  except  for  those  presented  in  [2]. 

Positive  definite  Toeplitz  matrices  play  an  important  role  in  several  aspects 
of  classical  analysis  including  the  classical  theory  of  Szego  polynomials,  the  tri¬ 
gonometric  moment  problem,  and  the  Caratheodory  coefficient  problem  [l,  16]. 
Several  Toeplitz  solvers  are  direct  manifestations  of  algorithms  and  formulas  that 
were  developed  in  the  classical  theory  (see,  e.g.,  [20,  19,  3]).  In  [3],  the  superfast 
Toeplitz  solver  presented  in  [17]  and  [22]  is  described  in  terms  of  a  generalization 
of  Sciiur’s  classical  algorithm. 

In  this  paper  we  give  a  general  description  of  the  Generalized  Schur  Algo¬ 
rithm  and  its  use  in  the  superfast  solution  of  Toeplitz  systems.  We  then  present 
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some  results  of  our  experimental  code  for  the  superfast  solution  of  real  positive 
definite  Toeplitz  matrices.  Refer  to  [3]  for  a  development  of  the  Generalized  Schur 
Algorithm  and  its  relation  with  fast  Toeplitz  solvers,  and  see  (2,4]  for  descriptions 
of  the  implementation  of  the  algorithm  for  complex  and  real  positive  definite  Toe¬ 
plitz  matrices. 


2.  Some  Toeplitz  Solvers 

The  Toeplitz  solvere  we  consider  can  be  divided  into  two  phases:  a  factorization 
phase  and  a  solution  phase.  In  this  section  we  summarize  some  fast  Toeplitz 
solvers  based  on  the  Szego  recursions  and  Schur’s  algorithm.  We  then  present  the 
generalization  of  Schur’s  algorithm  for  the  superfast  solution  of  a  positive  definite 
Toeplitz  system  of  equations. 

2.1.  The  Szego  Recursions.  The  Hermitian  positive  definite  Toeplitz  matrix 
M—  A/n+1  =  k=,  o  =  M*  defines  an  inner  product  on  the  set  of  polynomi¬ 

als  of  degree  at  most  n  by 

<\j ,\k  >  :  =  ^_* 

The  monic  polynomials  \j ,  deg\';  =j,  j  =  0,  ■  •  •  ,n  that  are  orthogonal  with 
respect  to  this  inner  product  are  the  (monic)  Szego  polynomials  determined  by  M . 
These  polynomials  can  be  constructed  using  roughly  n2  multiplications  and  addi¬ 
tions  by  the  Szego  recurrence  relations. 
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Szego  Recursions: 

Xq  :  =  1,  <50  •  =  ''o> 

for  .7  =  0,  1,  ...  ,  72—1 

lj+1  :  =  ~<1  ,\Xj>/&j 
Xy+i(X)  :  =  Xxy(X)  +  Ty+i  Xy(X) 

^y+i  :  =  <5y  (1-hy+i  P) 


where  Xy(X)  :  =  x(l/\)  is  the  polynomial  obtained  by  conjugating  and  revers¬ 

ing  the  coefficients  of  Xy-  The  quantities  7-  are  called  Schur  parameters ;  they 
serve  as  intermediate  quantities  in  the  Szego  recursions  and  are  often  of  physical 
and  mathematical  significance. 

The  Szego  recursions  comprise  the  first  phase  of  several  Toeplitz  solvers.  In 
fact,  they  are  equivalent  with  the  Levinson-Durbin  for  the  solution  of  the  Yule- 
Walker  equations  of  linear  prediction.  Define  R  to  be  the  unit  right  triangular 
matrix  whose  k th  column  contains  the  coefficients  of  x*  (0  <  k  <n).  By 
definition  of  the  Szego  polynomials  we  have 

R*MR=D,  M~l  =  R  D~l  R~* , 

where  D  :=  diag[^]0”  and  8k  :=  <Xk’Xk>-  Thus,  the  Szego  polynomials  deter¬ 
mine  the  reverse  Choleski  factorization  of  A/-1.  This  is  the  factorization  that 
underlies  the  solution  phase  of  Levinson’s  algorithm  for  the  solution  of  a  Toeplitz 
system  of  equations. 

In  Trench’s  algorithm  a  relationship  between  the  entries  of  M~l  is  used  to 
construct  M~l  from  its  last  column,  which  is  given  by  xn  an<^  8n.  This  relation¬ 
ship  follows  from  the  Christoffel-Darboux  formula  for  Szego  polynomials  [24,  20]. 
The  matrix  interpretation  of  the  Christoffel-Darboux-Szego  formula  is  the 
Gohberg-Semencul  formula 
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<5*  Mnh  =  ?’i  T*  -  To  To 

where 

0  0  •  •  •  0 

Po  0  •  .  ; 

:  o 

Pn- 1  •  •  '  Po  0 

1  Pn- 1  '  Po 

0  1 

Pn  —  1 

o  •  •  •  o  1 

and  Xn(^)  =  YjoPj  >  (P„  =  l)-  See,  e.g.,  [12].  Thus,  the  Gohberg-Semencul  for¬ 
mula  provides  a  factorization  of  M_1  that  is  determined  by  the  last  Szego  polyno¬ 
mial  Xn  and  its  squared  norm  8n.  The  product  M~lb  can  then  be  obtained  in 
0(n  log?r)  operations  by  using  Fast  Fourier  Transform  (FFT)  techniques.  This 
provides  an  0(n\ogn)  procedure  for  the  second  phase  of  a  Toeplitz  solver. 


2.2.  Schur’s  Algorithm.  Another  class  of  Toeplitz  solvers  is  based  on  the  clas¬ 
sical  algorithm  of  Schur  [23].  Let  <f  be  a  Schur  function ,  which  is  a  holomorphic 
mapping  of  the  open  unit  disk  in  the  complex  plane  into  its  closure.  Schur’s  algo¬ 
rithm  generates  a  continued  fraction  representation  of  a  given  Schur  function 
4>  =  (p0  by  the  successive  application  of  linear  fractional  transformations  (LFTs). 
Given  a  Schur  function  define  7n  =  <fin_ j(0)  and 


^n-l(X)  ~  In 
1  -Tn0n-l(X) 


It  is  easy  to  see  that  <fn  is  also  a  Schur  function.  Furthermore,  if  |7„  |  =  1  then 
(X)  =7n ;  in  this  case  0  is  a  rational  function  and  the  algorithm  terminates. 
Otherwise,  we  repeat  the  procedure  to  obtain  (fn  +  y  from  <f>n.  In  this  way  we 
obtain  the  Schur  continued  fraction  representation  of  , 
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<t>  =  7\«(<r>n)  =  ^lp<2°  ■  '  ■  °tn{<l>n)i 


where 


tn(T)  = 


In  +Xr 

1  +  7„  X  T 


The  function  Tn( 0)  is  referred  to  as  the  nt/i  approximant  of  (j),  and  <j)n  is  called  the 
nth  tail  of  4>.  Thus,  Schur’s  algorithm  effects  a  sequence  of  elementary  linear  frac¬ 
tional  transformations  on  the  initial  Schur  function  4>,  and  the  Schur  parameters 
•  are  intermediate  quantities  that  determine  the  elementary  linear  fractional 
transformations. 

To  implement  Schur’s  algorithm  we  write  each  Schur  function  as  the  quotient 

ft„  j 

of  formal  powrer  series,  <ftn  =  73 —  = - — ,  with  f3n  0  >  0  as  a  partial  nor- 

Ef-cA,  ,j 

realization.  Organizing  Schur’s  algorithm  so  that  the  coefficient  pairs  a0  k  ,80  k  are 
entered  and  processed  sequentially  leads  to  the  progressive  form  of  Schur  algo¬ 
rithm,  which  is  given  below. 


Schur’s  Algorithm: 

for  Ic  =  1,2,3,  •  while  | ')k  |  <  1 
enter  aQk_x,  (30k_l 
for  j  =  1,  2,  3,  •  •  •  ,  k- 1 


a 


i  ,k  —  j  —  1 


T*  =  a*-i,o/A-i,0’ 

0k. o  =  0k- 1,0  (1- hk  I2) 


L  "E 

-7;  1 


aj-\,k-j 

0j-\,k-j  J 


If  the  polynomials  a^n\  /3<jE  of  degree  less  than  n  formed  by  the  first  n  terms  of 
Q'0,  0q  are  input,  then  the  progressive  Schur  algorithm  calculates  Qjn~k\  l3[n~kl 
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and  ')k,  Ic  =  1,  ■  •  ■  ,  ?!  using  roughly  n2  multiplications  and  additions. 

In  contrast  with  the  Szego  recursions,  Schur’s  algorithm  provides  a  class  of 
fast  Toeplitz  solvers  based  on  the  Choleski  factorization  of  A 1 ,  M  =  L  D  L  * ,  with 
L  unit  lower  triangular  (see,  for  example,  [19,  22,  3]).  In  particular,  we  have  the 
following. 

Proposition  2.1.  If  M  =  [l^j-k}?  k- o  =  M*  is  positive  definite,  then 
Ocq  ■  =  —  /T  +  1,  /30j  =  Ji-,  j  =0,1,  .  .  .  ,n—  1,  are  the  first  n  coefficients  of 
formal  power  series  a0,(30  such  that  4>0  =  ct0/  is  a  Schur  function. 

Proposition  2.2.  Let  M ,  q^" \  be  as  in  Proposition  2.1.  Then 

1.  The  Schur  parameters  obtained  by  Schur’s  algorithm  applied  to  a$J \ 

are  identical  with  the  Schur  parameters  generated  by  the  Szego  re¬ 
cursions  applied  to  M . 

2.  The  lower  triangular  matrix  LD  =  [r;-  k\*,k-o  1S  given  by 
Tj,k  ~  Pk  ,j  —  k  —:j)- 

We  remark  that  the  connection  between  Schur  functions  and  Toeplitz  matrices  is 
provided  by  the  solution  of  the  Caratheodory  coefficient  problem  and  the 
correspondence  between  Caratheodory  functions  and  Schur  functions.  See,  for 
example.  [16, l].  Proposition  2.2  follows  from  by  comparison  of  the  resulting 
Szego  and  Schur  recursions. 

Schur’s  algorithm  can  also  be  formulated  in  terms  of  the  LFTs  Tn .  It  is 
shown  in  [2,  3]  that 


T»{r) 


Zn  +  Vn  T 

V  n  +  In  T  ' 


where  ,  t]n  are  polynomials  that  satisfy  the  recurrence  relations 


Vn-  \ 

£n-l 

In 

£o 

0 

_£n  —  1 

7?n—  1. 

.!  . 

1 

. 

and  £n(X):=Xn  £n(l/X).  r]n  (X)  :  =  X"  7?n  ( 1  /X).  We  refer  to  and  i]n  ,  which 
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have  degree  less  than  n  for  ?i  >  1,  as  the  nth  Sclivr  polynomials  associated  with 
the  Schur  function  4>0.  The  Schur  polynomials  are  generalizations  of  the  Schur 
parameters  in  the  sense  that  they  determine  compositions  of  the  elementary  linear 
fractional  transformations  that  the  Schur  parameters  determine. 

The  Schur  polynomials  are  also  important  in  that  they  can  be  used  to  obtain 
the  Szego  polynomials. 

Proposition  2.3.  Let  M,  fib’1'1  be  as  in  Proposition  2.1,  and  let  ik  ,t]k  be 

the  /t th  Schur  polynomials  determined  by  o^n),  /3^n\  Then  the  monic  Szego 
polynomials  corresponding  with  M  are  determined  by 

x*(X)  =  Vk  +  /  x. 


2.3.  The  Generalized  Schur  Algorithm.  Let 

Tn,k  =  +  1  0  +2°  ■  0  tn  +k 

be  the  LPT  generated  by  k  steps  of  Schur’s  algorithm  applied  to  Qn .  and  let 
k’Vn  ,k  be  Schur  polynomials  corresponding  with  Tn  k .  Note  that  for  any 
k  >  0  we  have  6n  =  Tn  k  (<f>n  +k ):  that  is,  the  k  th  tail  of  <t>n  is  the  {n+k)  th  tail 
of  6. 

The  Generalized  Schur  Algorithm  is  based  on  a  doubling  procedure  for  gen¬ 
erating  T  ..  fr< .m  7'.,.  It  can  be  described  in  general  terms  as  follows. 
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Generalized  Schur  Algorithm. 

To  obtain  £2h  .ri2n  from  Q^~n  Z, 

0.  Obtain  fn,  ??„  from  (3^n\ 

1.  Obtain  a{n\  /?£n)  from  a^n\  0^2n\  f  ,  rjn  using  4>n  =  T~\<bQ), 

2.  Obtain  0n  n,Vn  n  from  a{n\  0^  as  0n,Vn  were  obtained  from 

3.  Obtain  02n,V2n  from  fn>  *?«>  n^n.n  usine  r2n  =  ^0, „  “T'n.n- 


The  algorithm  is  started  by  performing  Step  0  directly  by,  for  example,  set¬ 
ting  =  cvoo/Z^oo  an<^  ??i  —  f-  More  generally,  we  can  use  Schur’s  algorithm  to 

generate  [7^]"°  for  a  small  value  of  n0,  and  obtain  £no.  t?no  from  the  Schur  param¬ 
eters  and  the  recursions  (2.1).  Now  suppose  Step  0  has  been  completed  for  some 
?!  >?;0. 

We  have 


n 


Qf 


T.-'M 


Q0Vn  ~  Po 0n 
00 Vn  ^0011 


and  in  [2.3]  it  is  shown  that 


Qo%  -  Potn  =  T„  +  ]  X"  +  0(X"  +  1), 


@oVn  ~^0n  =^nX"  +  0(M+1). 


We  can  therefore  take 

=  (a-oVn  ~  0O0n  )/  X"  ,  Pn  =  (/V?n  “  Qbfn  )/  X" 

These  equations  allow  us  to  obtain  ajM,  from  alfn \  Vn  in  SteP  L 

Step  2  is  the  doubling  step.  We  obtain  £n  n  and  ??n  „  from  a-Zn)  and  as 
0o  n  =  0n  a°d  Vo  n  —  Vn  Ave,'e  obtained  from  Qq11'1  and  B^n\ 
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Step  3.  the  composition  of  Tn  and  Tn  n ,  is  given  by 

^0,2 n  Wo,n  ,n  ^0,n  Vn  .n  • 

Wo, 2 n  ?0 ,n  ^>n  ,n  “b  Vo,n  Wn  .n  ■ 

Thus,  the  Generalized  Schur  Algorithm  consists  of  various  polynomial  multi¬ 
plications  and  additions  performed  in  a  recursive  manner.  The  multiplication  of 
polynomials  can  be  efficiently  performed  using  standard  FFT  techniques.  This 
results  in  an  0(nlog2n)  algorithm  for  computing  and  r\n ,  where  n  =  2".  The 
Toeplitz  system  of  equations  Mn  +  1  x  =  b  can  then  be  solved  by  forming  XT  from 
,  77 n  by  Proposition  2.3,  and  using  the  Gohberg-Semencul  formula  to  perform 
the  solution  phase  of  the  algorithm  in  0(nlog?i)  arithmetic  operations. 

The  recursive  nature  of  the  algorithm  is  easily  handled  in  a  programming 
language  that  allows  self-referencing  subroutines  (e.g..  PASCAL  or  APL).  In  such 
an  environment  it  is  easy  to  obtain  a  rudimentary  implementation  of  the  algo¬ 
rithm.  However,  in  order  to  obtain  a  FORTRAN  implementation  of  the  algo¬ 
rithm  we  must  consider  the  doubling  procedure  and  resulting  storage  requirements 
more  carefully. 

To  contrast  the  Generalized  Schur  Algorithm  with  Schur’s  algorithm,  the 
coefficients  of  ak,  ,3k  that  are  calculated  by  each  algorithm  are  displayed  in  Figure 
1  (for  n o  =  1 ).  While  the  classical  Schur  algorithm  generates  a[.n~k\  j3[n~k\  . 

0  <  k  <  n  from  Q’^n\  \  the  number  of  coefficients  of  ak,  (3k  calculated  by  the 

Generalized  Schur  Algorithm  is  equal  to  the  largest  power  of  two  that  divides  k . 
unless  this  power  of  two  is  less  than  n0.  Thus,  the  Generalized  Schur  Algorithm 
computes  only  pieces  of  the  Choleski  factorization  of  the  Toeplitz  matrix. 
Nevertheless,  we  can  obtain  every  Schur  parameter;  these  parameters  are  often  of 
physical  and  mathematical  significance. 


FIGURE  i 


The  A'th  row  of  squares  represents  the  coefficients  of  ak  and  0k  that  are  calculated  by  Schur’s  algo¬ 
rithm  and  the  Generalized  Schur  Algorithm  ( k  =0,...,n—  l).  The  filled  squares  represent  the  quanti¬ 
ties  that  are  calculated  during  one  iteration  of  each  algorithm.  In  Schur’s  algorithm  the 
coefficients  of  o0)  0o  are  processed  one  at  a  time.  In  the  Generalized  Schur  Algorithm  the 
coefficients  of  q0  ,/30  are  processed  in  segments,  each  segment  being  twice  the  size  of  the  previous 
segment. 
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3.  Numerical  Experiments 

We  have  a  working  FORTRAN  code  for  the  Generalized  Schur  Algorithm  for 
real  data  based  on  the  procedure  described  in  [4].  Our  code  implements  the  Gen¬ 
eralized  Schur  Algorithm  in  real  arithmetic,  and  uses  dual  split-radix  FFT  algo¬ 
rithms  to  avoid  the  explicit  use  of  the  bit  reversal  permutation.  The  code  uses 
6  n  4-  O(l)  real  storage  locations,  and  an  additional  3n/2  real  storage  locations 
for  the  Toeplitz  matrix  and  the  roots  of  unity  (sines  and  cosines)  required  by  the 
FFTs.  The  ?ith  monic  Szego  polynomial  xn  an<^  its  squared  norm  6n ,  where 
n  =2",  are  obtained  using  (8/3)n  log|n  +  0(nlog2n)  real  multiplications  and 
(l6/3)?i  log2?i  -f  0(nlog2n)  real  additions.  We  plan  to  make  our  code  available 
to  the  scientific  community  soon. 

In  [4]  it  is  shown  that  xn  <  K  can  be  obtained  using  less  than  8  n  log  |n  total 
real  arithmetic  operations.  To  our  knowledge,  this  is  the  smallest  operation  count 
obtained  for  a  superfast  Toeplitz  solver.  Since  the  Szego  recursions  require  more 
than  2  n2  operations,  the  Generalized  Schur  Algorithm  requires  less  computational 
work  for  n  =  2"  >  256.  However,  our  code  does  not  incorporate  the  considera¬ 
tions  discussed  in  [4,  Section  5].  The  operation  count  for  our  code  is  therefore 
somewhat  more  than  the  operation  count  of  [4]  in  the  lower  order  terms. 

When  the  Gaussian  or  Choleski  factorization  is  used  to  solve  the  system  of 
equations  Ax  —  b,  the  computed  solution  x  satisfies  (A+E)x  =  b  with 
l|£ II  <<<«’>)  IU  II,  where  e  is  the  machine  precision  and  <f>(n )  is  a  slowly  growing 
function  of  n.  It  follows  that  the  residual  r=b—Ax=Ex  and  the  error 
e  =  A~lE  x  satisfy 

llr  II  <  IP  ||  p||,  II  e  ||  <t<t>{n)  ||A  ||  ||A-1||  p||. 

Cybenko  [9,  10]  has  shown  that  4>{n)  is  of  order  n 2  for  the  Szego  recursions  (i.e., 
for  the  Levinson-Durbin  algorithm  for  the  Yule-Walker  equations).  As  the 
Choleski  decomposition  method  also  yields  4>{n)  of  order  n~,  we  can  conclude  that 
the  accuracy  of  the  Szego  recursions  is  comparable  with  that  of  Choleski’s 


-  12  - 


method.  Accurate  results  can  therefore  be  expected  from  the  Szego  recursions 
when  the  positive  definite  Toeplitz  matrix  is  well  conditioned. 

We  now  present  some  numerical  experiments  to  compare  the  performance  of 
the  Generalized  Schur  Algorithm  with  the  Szego  recursions.  We  are  pleased  to 
report  that  the  results  we  have  obtained  are  so  far  favorable.  In  particular,  the 
Generalized  Schur  Algorithm  displays  roughly  the  same  stability  characteristics 
(e.g.,  the  rates  of  growth  in  computed  residuals  and  errors)  as  the  Szego  recur¬ 
sions.  The  experiments  were  performed  on  the  VAX  11/750  at  Northern  Illinois 
University. 

We  performed  our  experiments  using  the  leading  principal  submatrices  of  two 
classes  of  real  positive  definite  Toeplitz  matrices  A/n+1  =  [/4;_*]y  k  _ 0  =  . 

We  refer  to  these  classes  as  type  A  and  type  B. 

*2 

A(n,9)  refers  to  A/n  +  1  where  p,j  —  9}  ,  — 1<#<1.  These  matrices  have  a  par¬ 
ticularly  special  structure  [14].  The  Schur  parameters  are  given  by  7  -  =  (—#)+ 
For  fixed  n  these  matrices  become  ill-conditioned  as  9  tends  to  ±1.  For  fixed  9 
the  condition  numbers  tend  toward  a  limiting  value  as  n  increases.  These  are 
numerically  banded  matrices,  where  the  bandwidth  depends  on  6  and  the  smallest 
positive  machine  number. 

B{n  ,9x,92,p)  refers  to  the  autocorrelation  matrix  Mn  +  1  of  order  n+1  of  a  vec¬ 
tor  x  =  [£;]oI6n_I  given  by 

n /  4 

=  vj  +  £  cos (m), 

*-1 

where  th  uik  are  n/4  equispaced  integral  multiples  of  7r/ 8 n  in  [^,^2] 
(0<*i<*2<"/2),  and  rjj  is  randomly  generated  in  [— p,p\.  In  other  words, 
A/n+1  =  Tt  T,  where  T  is  the  Toeplitz  matrix  formed  by  the  first  n+1  columns 
of  the  circulant  matrix  whose  first  column  is  x.  These  matrices  become  ill- 
conditioned  as  p  tends  to  zero,  and  the  ill-conditioning  accelerates  as  62  —  9 , 
decreases. 


We  first  compare  the  accuracy  of  the  Szego  recursions  (Sz)  and  the  General¬ 
ized  Schur  Algorithm  (GS)  for  computing  the  monic  Szego  polynomial  xn  and  the 
Schur  parameters  [7 -j".  In  Tables  1  and  2  we  give  the  square  root  of  the  sum  of 
the  squared  errors  (the  2-norm  of  the  error)  as  well  as  the  maximum  absolute 
error  in  the  coefficients  of  xn  and  [lj } "  for  the  indicated  matrices  of  type  A  and 
type  B.  The  results  of  the  Szego  recursions  in  double-precision  arithmetic  are  used 
as  exact  answers  for  the  error  calculation.  Note  that  the  errors  for  the  two  algo¬ 
rithms  seem  to  grow  at  roughly  the  same  rate  as  n  increases.  In  well-conditioned 
problems,  the  Generalized  Schur  Algorithm  typically  results  in  a  somewhat  larger 
error  than  the  Szego  recursions,  but  usually  this  difference  is  less  than  an  order  of 
magnitude.  On  the  other  hand,  the  Generalized  Schur  Algorithm  occasionally 
yielded  more  accurate  results  in  the  ill-conditioned  problems  of  type  B. 


TABLE  1 
Phase  One  Errors 


A  (4096,0  6) 

Szego  polynomial  errors  Schur  parameter  errors 

n  2-norm _  max  2-norm  max 


GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0  184e-06 

0.219e-06 

0.119e-06 

0.149e-06 

0  802e-07 

0.112e-06 

0  596e-07 

0  894e-07 

8 

0  315e-06 

0  635e-06 

0.209e-06 

0.328e-06 

0  136e-06 

0  213e-06 

0.782e-07 

0  104e-06 

16 

0.353e  06 

0.738e-06 

0.201e-06 

0.328e-06 

0  181e-06 

0  227e-06 

0.782e-07 

0  104  e- 06 

32 

0  557e-06 

0  780e-06 

0.209e-06 

0  343e-06 

0.229e-06 

0.227e-06 

0  785e-07 

0  104  e- 06 

64 

0  719e-06 

0  780e-06 

0.179e-06 

0.343e-06 

0  297^06 

0.227e-06 

0.988e-07 

0  104e-06 

128 

0  983e-06 

0  780e-06 

0  184e-06 

0  343e-06 

0  371e-06 

0.227e-06 

0  988e-07 

0. 104  e- 06 

256 

0.120e-05 

0  780e-06 

0.238e-06 

0.343e-06 

0.436e-06 

0.227e-06 

0.988e-  07 

0  104e-06 

512 

0.149e-05 

0.780e-06 

0.358e-06 

0.343e-06 

0.512e-06 

0.227e-06 

0.988e-07 

0  104e-06 

1024 

0. 175e-05 

0.780e-06 

0.477e-06 

0.343e-06 

0.579e-06 

0.227e-06 

0.988e-07 

0  104e-06 

2048 

0.199e-05 

0.780e-06 

0.596e-06 

0.343e-06 

0.639e-06 

0.227e-06 

0  988e-07 

0  104e-06 

4096 

0.225e-05 

0.780e-06 

0.775e-06 

0.343e-06 

0.698e-06 

0.227e-06 

A  AO O  -  A-* 
U.JOOL'Jl 

0  104  e- 06 

n 

2-: 

A  (4096,0  75) 

Szego  polynomial  errors 

norm  max  2-] 

Schur  parameters  errors 
norm  max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.333e-06 

0.720e-06 

0.238e-06 

0.536e-06 

0.942e-07 

0.321e-06 

0.894e-07 

0.298e-06 

8 

0  300e-05 

0  112e-04 

0.155e-05 

0  572e-05 

0  780e-06 

0.168e-05 

0  745e-06 

0.991e-06 

16 

0.475e-04 

0  276e-04 

0.187e-04 

0.111e-04 

0.473e-05 

0.249e-05 

0.205e-05 

0.991e-06 

32 

0  521e-04 

0.292e-04 

0.197e-04 

0.113e-04 

0.502e-05 

0.252e-05 

0.205e-05 

0.991e-  06 

64 

0.591e-04 

0  293e-04 

0.197e-04 

0  113e-04 

0.6l4e-05 

0.252e-05 

0.205e-05 

0  991e-06 

128 

0  696e-04 

0  293e-04 

0.197e-04 

0.113e-04 

0  736e-05 

0.2S2e-05 

0.205e-05 

0  991  e- 06 

256 

0  835e-04 

0.293e-04 

0.197e-04 

0  113e-04 

0  866e-05 

0.252e-05 

0  205e-05 

0  991e-06 

512 

0  986e-04 

0  293e-04 

0.198e-04 

0.1 13e-04 

0  100e-04 

0  252e-05 

0  205e-05 

0  991e-06 

1024 

0  120e-03 

0.293e-04 

0  198e-04 

0.113e-04 

0  118e-04 

0.252e-05 

0  205e-05 

0  991e-06 

2048 

0  140e-03 

0.293e-04 

0  198e-04 

0.1 13e-04 

0  135e-04 

0  252e-05 

0.205e-05 

0  991e-06 

4096 

0. 155e-03 

0.293e-04 

0  199e-04 

0.113e-04 

0  148e-04 

0  252e-05 

0  205e-0S 

0  991e-06 

n 

Szego  poly 
2- norm 

A  (4096,0. 

nomial  errors 

max 

83) 

2-; 

Schur  parameters  errors 
norm  max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0  449e-05 

0.190e-05 

0  286e-05 

0  143e-05 

0  136e-05 

0  786e-06 

0  131e-05 

0.775e-06 

8 

0  369e-04 

0  173e-03 

0.224e-04 

0  918e-04 

0.420e-05 

0.147e-04 

0.229e-05 

0  115e-04 

16 

0  180e-02 

0.290e-02 

0.785e-03 

0.114e-02 

0  876e-04 

0.775e-04 

0  643e-04 

0  321e-04 

32 

0.231e-01 

0.636e-02 

0  598e-02 

0.189e-02 

0  474e-03 

0.1 12e-03 

0  157e-03 

0.321e-04 

64 

0  309e-01 

0  655e-02 

0.678e-02 

0.190e-02 

0.638e-03 

0.1 13e-03 

0.159e-03 

0.321e-04 

128 

0  444e-01 

0  655e-02 

0  678e-02 

0.190e-02 

0  894e-03 

0  113e-03 

0  159e-03 

0  321e-04 

256 

0  547e-01 

0  655e-02 

0  679e-02 

0.190e-02 

0.109e-02 

0.113e-03 

0.159e-03 

0  321e-04 

512 

0  629e-01 

0  655e-02 

0.680C-02 

0.190e-02 

0  122e-02 

0.113e-03 

0  159e-03 

0  321e-04 

1024 

0  682e-01 

0  655e-02 

0  681e-02 

0  190e-02 

0  131e-02 

0.113e-03 

0  159e-03 

0  321e-04 

2048 

0  725e-01 

0.655e-02 

0  681e-02 

0  190e-02 

0  140e-02 

0  113e-03 

0.159e-03 

0  321e-04 

4096 

0  769e-01 

0  655e-02 

0  682e-02 

0.190e-02 

0.148e-02 

0.113e-03 

0  159e-03 

0  321e-04 

TABLE  2 
Phase  One  Errors 


£(4096,20  ‘,160  ',10  0) 

Szego  polynomial  errors  Schur  parameters  errors 

n  2-norm _  _ max _  2-norm  max 


GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0  507e-08 

0  314e-07 

0.301e-08 

0  310e-07 

0.999e-08 

0  999e-08 

0  952e-08 

0.952e-08 

8 

0  380e-07 

0.149e-06 

0.659e-07 

0.838e-07 

0  357e-07 

0  933e-07 

0  242e-07 

0  838 e- 07 

16 

C.422e-06 

0.850e-06 

0  185e-06 

0.411e-06 

0  134e-06 

0  287e-06 

0  888e-07 

0  I58e-06 

32 

0  776e-06 

0  878e-06 

0.294e-06 

0.404e-06 

0.208e-06 

0  340e-06 

0.888e-07 

0  158e-06 

64 

0  186e-05 

0.918e-06 

0.532e-06 

0.387e-06 

0  499e-06 

0  413e-06 

0  151  e- 06 

0  158e-06 

128 

0.256e-05 

0.111e-05 

0  841e-06 

0.504e-0o 

0  587e-06 

0  434e-06 

0  151e-06 

0  158e-06 

256 

0  399e-05 

0  138e-05 

0  138e-05 

0.692e-06 

0.683e-06 

0  477e-06 

0  151e-06 

0  158e-06 

512 

0  670e-05 

0. 181e-05 

0.238e-05 

0.734e-06 

0  781e-06 

0  540e-06 

0  151e-06 

0.158e-06 

1024 

0.11Se-04 

0  246e-05 

0  441e-05 

0.955e-06 

0  104e-05 

0  656e-06 

0  151e-06 

0  158e-06 

204P 

0.205e-04 

0  344e-05 

0.81  le-05 

0.966e-06 

0.139e-05 

0,919e-06 

0.151e-06 

0  158e-06 

4095 

0.854e-04 

0.270e-04 

0.297e-04 

0.1 19e-04 

0  547e-04 

0.168e-04 

O.250e-04 

0.616e-05 

£(4096,20  ',160  *,0.5) 

Szego  polynomial  errors  Schur  parameters  errors 


n 

2-norm 

max 

norm 

max 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.346e-07 

0.S87e-07 

0.257e-07 

0.515e-07 

0.328e-07 

0.570e-07 

0  239e-07 

0  515e-07 

8 

0  945e-07 

0.474e-06 

0  596e-07 

0.232e-06 

0.625e-07 

0  261e-06 

0.527e-07 

0  232e-06 

16 

0.807  e-05 

0.377e-04 

0.284e-05 

0  144e-04 

0.1 16e-05 

0  527e-05 

0  873e-06 

0  430e-05 

32 

0  llle-02 

0  480e-02 

0  406e-03 

0  195e-02 

0.467e-04 

0.205e-03 

0  204e-04 

0.937e-04 

64 

0.436e-02 

0  499e-02 

0  172e-02 

0.203e-02 

0.148e-03 

0  312e-03 

0  588e-04 

0  999e-04 

128 

0  501e-02 

0  524e-02 

0.156e-02 

0  187e-02 

0.185e-03 

0  362e-03 

0  588e-04 

0  999e-04 

256 

0  666e-02 

0.848e-02 

0.153e-02 

0.163e-02 

0  238e-03 

0.490e-03 

0.588e-04 

0  9 99 e- 04 

512 

0  796e-02 

0.182e-01 

0  158e-02 

0.170e-02 

0.280e-03 

0  809e-03 

0  588e-04 

0  999 e- 04 

1024 

0  139e-01 

0.316e-01 

0  189e-02 

0.231e-02 

0  428e-03 

0  140e-02 

0.588e-04 

0  999e-04 

2048 

0  461e-01 

0  659e-01 

0  303e-02 

0.353e-02 

0.145e-02 

0.311e-02 

0  116e-03 

0  173e-03 

4096 

0  506e-01 

0  413e+00 

0  131e-01 

0.104e+00 

0  376e-01 

0  221e+00 

0  132e-01 

0  830e-01 

£(4096,50  ',100  ',10) 

Szego  polynomial  errors  Schur  parameters  errors 

n _  _ 2-norm _  _ max  2-norm  max 


GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.393e-05 

0.409e-05 

0  238e-05 

0.262e-05 

0.169e-05 

0  177e-05 

0  161e-05 

0.173e-05 

8 

0  490e-04 

0  393e-04 

0.312e-04 

0.277e-04 

0  270e-04 

0.255e-04 

0  206e-04 

0  187e-04 

16 

0.124e-03 

0.307e-03 

0.830e-04 

0.154e-03 

0.585e-04 

0  120e-03 

0.427e-04 

0  631e-04 

32 

0  319e-03 

0  500e-03 

0.140e-03 

0.170e-03 

0  129e-03 

0.224e-03 

0.480e-04 

0  757 e- 04 

64 

0  411e-03 

0.647e-03 

0  120e-03 

0.165e-03 

0.177e-03 

0.333e-03 

0  480e-  04 

0  757 e- 04 

128 

0  644e-03 

0  890e-03 

0  141e-03 

0.197e-03 

0  245e-03 

0  479e-03 

0.634e-04 

0  757 e- 04 

256 

0  952e-03 

0  129e-02 

0  165e-03 

0  246e-03 

0  346e-03 

0  676e-03 

0  634e-04 

0  783 e- 04 

512 

0  152e-02 

0.231e-02 

0  196e-03 

0.309e-03 

0,530e-03 

0  104e-02 

0  634e-04 

0  989e-04 

1024 

0  260e-02 

0  325e-02 

0.231e-03 

0  357e-03 

0  904  e- 03 

0.155e-02 

0  726e-04 

0  114e-03 

2048 

0  303e-02 

0  439e-02 

0  254e-03 

0.396e-03 

0  122e-02 

0,217e-02 

0.726e-04 

0  114e-03 

4096 

0  395e-02 

0  641e-02 

0  282e-03 

0  436e-03 

0.173e-02 

0  312e-02 

0  769e-04 

0  114e-03 

-  14  - 


Our  second  accuracy  experiment  compares  the  Szego  recursions  and  the  Gen¬ 
eralized  Schur  Algorithm  for  the  solution  of  a  Toeplitz  system  of  equations.  The 
Gohberg-Semencul  formula  is  used  for  the  solution  phase  for  both  algorithms. 

Let  || A/ 1| j  and  ||AL||2  denote  the  norm  of  the  matrix  M  that  is  induced  by 
the  vector  1-norm  and  vector  2-norm,  respectively.  Then 

max  ||A^  e,  ||2  <  ||M||2  <  Vn"  UmHj.  (3.1) 

1  <  j  <  n 

Let  M  —  Mn ,  and  define  R  :  =  R  D~'k  so  that  M~x  =  R  R* .  It  is  shown  in  [2] 
that 

<I|M-'||2<..  ||«ej?-  (3.2) 

This  gives  bounds  on  the  2-norm  of  M~l  in  terms  of  the  last  column  of  R  (i.e.,  in 
terms  of  the  coefficients  of  the  normalized  Szego  polynomial  Xn-i/V^n-i  )• 

After  computing  a  solution  x  of  the  positive  definite  Toeplitz  system,  the 
residual  r  =  b  —  M  x  is  calculated  in  double  precision,  and  the  system  M  e  =  r 
is  solved  (in  0(n  lg  n )  operations)  to  obtain  an  estimate  e  of  the  error.  Finally, 

the  upper  and  lower  bounds  for  the  two  ratios  ipr  —  — — ^-|| — —  and 

||M||||f|| 

ip*  ~  ~rr — n~. j'—  7 .. — 7T  are  obtained  using  (3.1)  and  (3.2).  In  order  for  an  algo- 

m  \\m-'\\  iuii 

rithm  to  be  stable,  these  ratios  should  increase  slowly  with  n.  In  Tables  3  and  4 
we  display  the  computed  lower  bounds  on  the  computed  relative  residu¬ 

als  and  relative  errors,  and  upper  bounds  on  the  ratios  ipr  and  These  results 
were  obtained  using  6  =  [j  +l]/«o-  Similar  results  were  obtained  when  the  com¬ 
ponents  of  6  were  set  equal  to  one  and  when  they  were  randomly  generated.  The 
Szego  recursions  were  stopped  at  n  =  213. 

Note  that  the  sizes  of  the  residuals  and  errors  for  each  algorithm  are  roughly 
the  same,  as  are  the  rates  of  growth  in  these  quantities.  The  upper  bounds  on  i/v 
and  il>e  are  slowly  growing  for  both  algorithms.  The  errors  and  residuals  for  the 


TABLE  3 

Toeplitz  Solver  Residuals  and  Errors 


A  (21S,0.6) 


n 

IKA  11 2 

Hr 

Jb 

II 

H 

IU 

Hi 

II 

il _ 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.449?+01 

0.376?-06 

0.580e-06 

0.353?-06 

0.483?-06 

0.376?-06 

0.580?- 06 

0.594?- 07 

0.812?-07 

8 

0.493?+01 

0.221e-06 

0.615?-06 

0.2S3?-06 

0.631?-06 

0.266e-06 

0.740?- 06 

0.388?- 07 

0.967?- 07 

16 

0.495?+01 

0.383e-08 

0.444e-06 

0.594?-06 

0.898?-06 

0.54S?-06 

0.631?- 06 

0.906e-07 

0.137?- 06 

32 

0.495?+01 

0.280e-06 

0.457e-06 

0.540?-06 

0.67 7?- 06 

0.449?-06 

0.731?- 06 

0.825?- 07 

0.1 03?- Of. 

6-1 

0.49S?+01 

0.725e-06 

0.117e-06 

0.310e-06 

0.711e-06 

0.125?-05 

0.202?- 06 

0.139?- 06 

0.109?- 06 

128 

0.495?+01 

0.567e-06 

0.296e-06 

0.120?-05 

0.107?-05 

0.102e-05 

0.529?- 06 

0.184  ?-06 

0.163?- 06 

256 

0.495?+01 

0.105e-05 

0.205e-06 

0.171e-0S 

0.140?-05 

0.192e-05 

0.375?- 06 

0.261e-06 

0.214?- 06 

512 

0.495e+01 

0.150e-OS 

0.327e-06 

0.161e-05 

0.697e-06 

0.278?-05 

0.604e-06 

0.246e-06 

0.106?- 06 

1021 

0.495?+01 

0.206e-05 

0.399e-06 

0.222f-05 

0.893?-06 

0.383e-0S 

0,743?- 06 

0.338e-06 

0.136?- 06 

2018 

0.49Se+01 

0.288e-05 

0.446c-06 

0.299e-05 

0.934?-06 

0.538e-05 

0.833e-06 

0.457?- 06 

0.143?- 06 

4096 

0.495e+01 

0.410«-05 

0.501e-06 

0.418?-05 

0.944e-06 

0.767?-05 

0.937e-06 

0.638?- 06 

0.144?- 06 

8192 

0.495e+01 

0.595e-05 

0.683e-06 

0.600?-05 

0.109e-05 

0.111e-04 

0.128e-05 

0.916?- 06 

0.167?- 06 

16384 

0.495?+01 

0.571e-05 

O.S90?-OS 

0.107e-04 

0.901?- 06 

32768 

0.495e+01 

0.723e-05 

0.738e-05 

0.13S?-04 

0.113?- 05 

A  (216,0.75) 


IUC+1,  ||2 

Hr 

II 

IU 

II 

1 A 

,/• 

n 

Jb 

If 

II  i 

II 

Vf 

'  C 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.279?+02 

0.703?-06 

0.122?-0S 

0.434?-06 

G,lb9?-05 

0.321e-06 

0.557?- 06 

0.102?- 07 

0.397?-07 

8 

0.647?+02 

0.704?-06 

O.591e-0o 

0.278?-OS 

0.514?-05 

0.3l6?-06 

0.265?- 06 

0.281?- 07 

0.520?-07 

16 

0.799?-?02 

0. 170?-05 

0.163?-05 

0.905?-05 

0.115^04 

0.967?-06 

0.923?- 06 

0.741?- 07 

0.942?-  07 

32 

0.803?+02 

0.123?-05 

0.425?-05 

0.140?-04 

0.117?-04 

0.959?-06 

0.331?-0S 

0.114?-  06 

0.954?- 07 

64 

0.803?+02 

0.386e-05 

0.285?-05 

0.161?-04 

0.142?-04 

0.402?-0S 

0.296?-05 

0.131?- 06 

0.116?- 06 

128 

0.803?+02 

0.513?-05 

C.293?-05 

0.289?-04 

0.260?-04 

0.681?-05 

0.389?- 05 

0.235?- 06 

0.212?- 06 

256 

0.303?+02 

0.922?- 05 

0.2S5?-05 

0.52 0?-04 

0.504e-04 

0.148e-04 

0.408?- 05 

0.424?- 06 

0.411?- 06 

512 

0.803?+02 

0.384 ?-0S 

0.286?-05 

0.267?-04 

0.207 ?-04 

O.699?-0S 

0.521e-05 

0.217?- 06 

0.169?-  06 

1024 

0.803?+02 

0.610?-05 

0.300?- 05 

0.372?-04 

0.3 15?- 04 

0.120?-04 

0.590?-05 

0.303?- 06 

0.257?-  06 

2048 

0.803?+02 

0.891e-05 

0.288?- 05 

0.363?-04 

0.355?-04 

0.183?-04 

0.593?- 05 

0.296?- 06 

0.289?-  06 

4096 

0.803?+02 

0.951?-05 

0.322?-05 

0.317?-04 

0.303?-04 

0.200?-04 

0.679?- 05 

0.258?- 06 

0.246?- 06 

8192 

0.803?+02 

0.692?-05 

0.334?-05 

0.338?-04 

0.341?-04 

0.148?-04 

0.7I4?-05 

0.275?- 06 

0.278?- 06 

16384 

0.803?+02 

0.107?-04 

0.744?-04 

0.231?-04 

0.606?- 06 

32768 

0.803e+02 

0. 159?-04 

0.738?-04 

0.344e-04 

0.601?- 06 

A  (215,0.8) 


IlK ~J,  Hz 

Hr 

II 

IU 

II 

xh 

x h 

n 

lu 

II 

II  i 

II 

Vr 

V  t 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.721?+02 

0.161?-05 

0.26 5?- 05 

0.74’?-06 

0.1 29?-  05 

0.448?- 06 

0.739?- 06 

0.636?- 08 

0.1 10?- 07 

8 

0.320?+03 

0.461?-05 

0.546?-05 

0.661?-05 

0.213?-0S 

0.932?-06 

0.110?-05 

0.127?- 07 

0.408?- 08 

16 

0.610?+03 

0.462?-05 

0.124?-04 

0.327 ?-04 

0.826?-05 

0.985?-06 

0.265?-05 

0.329?- 07 

0.831?- 08 

32 

0.642?+03 

0.764?-05 

0.213?-04 

0.137?-03 

0.190?-04 

0.225?-05 

0.628?- 05 

0.131?- 06 

0.182?-07 

64 

0.642?+03 

0.112?-04 

0.192?-04 

0.162?-03 

0.400?-04 

0  4  68?- 05 

0.800?- 05 

0.155?- 06 

0.383?- 07 

128 

0.642?+03 

0.427?-05 

0.1S4?-04 

0.190?-03 

0.123?-03 

0.249?-05 

0.900?- 05 

0.181?- 06 

0.118?-  06 

256 

0.642?+03 

0.126?-04 

0.172?-04 

0.27 9?- 03 

0.291?-03 

0.101?-04 

0.138?-04 

0.267?- 06 

0.278?- 06 

512 

0.642?+03 

0. 113?-04 

0.174?-04 

0.264?-03 

0.129?-03 

0.121?-04 

0.186?- 04 

0.253?- 06 

0.123?- 06 

1024 

0  642?+03 

0.876?-05 

0.172?-04 

0.363?-03 

0.264?-03 

0.120?-04 

0.237?- 04 

0.347?- 06 

0.253?- 06 

2048 

0.642?+03 

0.486?-05 

0.175?-04 

0.392?- 03 

0.301?-03 

0.812e-05 

0.292?-04 

0.375?- 06 

0.288?- 06 

4096 

0.642?+03 

0.817e-05 

0.150?-04 

0.367 ?-03 

0.295?-03 

0.156?-04 

0.288?- 04 

0.351?- 06 

0.282?- 06 

8192 

0.642?+03 

0.116?-04 

0.178?-04 

0  421?-03 

0.372?-03 

0.241?-04 

0.371?- 04 

0.403?-  06 

0.356?- 06 

16384 

0.642?+03 

0.221?-04 

0.891?-03 

0.482?-04 

0.852?- 06 

327G8 

0.642?+03 

0.204?-04 

0.906?-03 

0.458?-04 

0.867?- 06 

TABLE  4 

Toeplitz  Solver  Residuals  and  Errors 


£(214,20  160  ’,50) 


n 

IIa/.-A  ll2 

l!r  11 

_ El _ 

- 

lull 

m  ii 

0. 

v/;, 

r 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.133e+01 

0.824*07 

0.824e-07 

0.760e-07 

0.760e-07 

0.545e-07 

0.545e-07 

0.538e-07 

0.538*07 

8 

0.178e+01 

0.211*06 

0.235*- 06 

0.136*- 06 

0.249*- 06 

0.828 *-07 

0.922*07 

0.701e-07 

0.129*- 06 

16 

0.192e+01 

0.372e-06 

0.221*- 06 

0.1 12e-06 

0.170e-06 

0.763e-07 

0.453e-C7 

0.529*-07 

0.806e-07 

32 

0.200e+01 

0.515*- 06 

0.305*06 

0.468*- 06 

0.173e-06 

0.889e-07 

0.526e-07 

0.212e-06 

O  '  7  *-07 

64 

0.203e+01 

0.408e-06 

0.248e-06 

0.198e-06 

0.878e-07 

0.657*- 07 

0.399*- 07 

0.883e-07 

C  91*07 

128 

0.205e+01 

0.538e-06 

0.389*  06 

0.270*-  06 

0.224e-06 

0.840e-07 

0.607*- 07 

0.119e-06 

0.988e-07 

256 

0.206e+01 

0. 101*- 05 

0.736*- 06 

0.904*- 06 

0.650e-06 

0.1S6e-06 

0.114e-06 

0.396*-06 

0.285*  06 

512 

0.206e+01 

0.395*- 05 

0.771e-06 

0.393*05 

0.677*- 06 

0.612*- 06 

0.120*- 06 

0.172e-05 

0.296*  06 

1024 

0.207e+01 

0.115*04 

0.117e-05 

0.115e-04 

0.110*05 

0.175e-05 

0.177e-06 

0.501e-05 

0.482e-06 

2048 

0.208e+01 

0.228 e-04 

0.514e-05 

0.229 e-04 

0.514*05 

0.339e-05 

0.763e-06 

0.992*05 

0.223e-05 

4096 

0.209e+01 

0.501*- 04 

0.164e-04 

0,501e-04 

0.164e-04 

0.723e-05 

0.237e-05 

0.216*04 

0.708*-  05 

8192 

0.212e+01 

0.102e-03 

0.375e-04 

0.102e-03 

0.375*- 04 

0.144*- 04 

0.532e-05 

0.432*-04 

0.160*- 04 

16384 

0.217e+01 

0.196e-03 

0.197e-03 

O.OOOe+OO 

0.265*- 04 

O.OOOe+OO 

0.814*04 

O.OOOe+OO 

32768 

0.614e+01 

0.469e-03 

0.470*03 

0.000e+00 

0.399*- 04 

0.000e+00 

0.520*-04 

O.OOOe+OO 

n 

lUCA  Hz 

llrll 

_ UH _ 

B(  21S 

,20  160  \5) 

IUII 

IUII 

0 

r 

0, 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.1S8e+01 

0.847*- 07 

0.942*- 07 

0.964*07 

0.119*06 

0.484*07 

0.539*07 

0.560*07 

0.689*  07 

8 

0.382e+01 

0.304e-06 

0.262*- 06 

0.192*06 

0.171*06 

0.732*07 

0.631*07 

0.448*07 

0.398*07 

16 

0.440e+02 

0.486*05 

0.314*05 

0.338*05 

0.559*06 

0.137*06 

0.886*  07 

0.676*07 

0.112*07 

32 

0.974e4-02 

0.273e-04 

0.272*04 

0.201*04 

0.112*04 

0.696*07 

0.693*  07 

0.180*06 

0.100*05 

64 

0.112e+03 

0.391*04 

0.285*  04 

0.266*04 

0.161*04 

0.810*07 

0.590*07 

0.207*06 

0.125*06 

128 

0.118e+03 

0.423e-04 

0.363*04 

0.299*04 

0.194*04 

0.800*  07 

0.687*07 

0.220*06 

0.143*06 

256 

0.121e+03 

0.600*  04 

0.426*04 

0.471*04 

0.201*04 

0.109*06 

0.777*07 

0.337*06 

0.144*06 

512 

0.123e+03 

0.729*04 

0.457*04 

0.662*04 

0.196*04 

0.132*06 

0.827*07 

0.468*06 

0.138*06 

1024 

0.124e+03 

0.914*04 

0.484*04 

0.830*04 

0.189*04 

0.161*06 

0.855*07 

0.582*06 

0.133*06 

2048 

0.124e+03 

0.933*04 

0.456*04 

0.833*  04 

0.144*04 

0.160*06 

0.783*07 

0.580*06 

0.101*06 

4096 

0.126e+03 

0.675*- 04 

0.472*04 

0.500*  04 

0.705*05 

0.113*06 

0.787*07 

0.345*06 

0.487*07 

8192 

0.128*4-03 

0.607e-04 

0.753*04 

0.235*  04 

0.530*04 

0.996*07 

0.124*06 

0.160*06 

0.360*06 

16384 

0.132e+03 

0.924e-04 

0.703*04 

0.000e+00 

0.144*06 

0.000t+00 

0.462*06 

O.OOOe+OO 

32768 

0.535e4-03 

0.922e-03 

0.129*03 

0.000e+00 

0.850*06 

O.OOOe+OO 

0.148*06 

O.OOOe+OO 

B(21S,70  110  •,  1.0) 


IU4"+i  lls 

IU 

II 

IU 

II 

0, 

" 

_ UM! _ 

_ 1U1L 

GS 

Sz 

GS 

Sz 

GS 

Sz 

GS 

Sz 

4 

0.276e+03 

0.223*04 

0.188*04 

0.513*05 

0.109*04 

0.717*07 

0.606*  07 

0.113*07 

0.239  *  07 

8 

0.771e+04 

0.654*03 

0.147*02 

0.273*03 

0.125*03 

0.482*07 

0.109*06 

0.183*07 

0.840*08 

16 

0.819e+04 

0.720*03 

0.148*02 

0.332*03 

0.238*03 

0.342*07 

0.704*07 

0.200*  07 

0.143*07 

32 

0.872e+04 

0.118*02 

0.161*02 

0,757*03 

0.662*03 

0.505*07 

0.687*07 

0.417*07 

0.364*07 

64 

0.926e+04 

0.222*02 

0.179*02 

0.141*02 

0.986*03 

0.905*07 

0.728*07 

0.722*07 

0.504*07 

128 

0.939e+04 

0.331*02 

0.208*02 

0.270*02 

0.124*02 

0.132*06 

0.831*07 

0.136*06 

0.619*07 

256 

0.953e+04 

0.419*02 

0.209*02 

0.394*02 

0.122*02 

0.165*06 

0.821*07 

0.194*06 

0.600*  07 

512 

0.959e+04 

0.584*02 

0.196*02 

0.554*02 

0.994*03 

0.226*06 

0.755*07 

0.271*06 

0.483*07 

1024 

0.963e+04 

0.689*02 

0.201*02 

0.645*02 

0.888*03 

0.254*06 

0.759*07 

0.313*06 

0.430*07 

2048 

0.969e+04 

0.749*02 

0.202*02 

0.720*02 

0.826*03 

0.283*06 

0.759*07 

0.348*06 

0.398*07 

4096 

0.978e+04 

0.825*02 

0.194*02 

0.785*02 

0.739*03 

0.318*06 

0.743*07 

0.375*06 

0.352*07 

8192 

0.101e+05 

0.831*02 

0.178*02 

0.117*01 

0.715*03 

0.323*06 

0.684*07 

0.544*06 

0.336*  07 

16384 

0.103e+05 

0.302*02 

0.789*02 

0.333*06 

0.358*06 

32768 

0.116e+05 

0.804*02 

0.277*01 

0.400*06 

0.111*05 
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Generalized  Schur  Algorithm  were  typically  larger  than  those  of  Szego  recursions 
in  the  well-conditioned  problems,  while  the  Generalized  Schur  Algorithm  some¬ 
times  performed  slightly  better  in  the  relatively  ill-conditioned  problems.  These 
differences  were  rarely  more  than  an  order  of  magnitude. 

Thus,  in  the  examples  we  investigated  the  accuracy  of  the  Generalized  Schur 
Algorithm  closely  followed  that  of  Szego  recursions.  Of  course,  while  we  did  not 
encounter  any  problems  where  the  accuracy  of  the  two  algorithms  was  greatly 
different,  more  experimentation  is  needed  to  obtain  a  better  confidence  in  the 
results  of  the  Generalized  Schur  Algorithm. 


Table  5  shows  average  CPU  times  required  by  the  Generalized  Schur  Algo¬ 
rithm  and  Szego  recursions  for  the  calculation  of  the  nth  Szego  polynomial  and 
the  n  Schur  parameters.  Note  that  the  Szego  recursions  would  have  used  about  36 
CPU  haul's  for  n=216,  compared  with  about  38  minutes  for  the  Generalized 
Schur  Algorithm.  The  dramatic  relative  efficiency  of  the  Generalized  Schur  Algo¬ 
rithm  for  large  problems  and  the  potential  reliability  of  the  algorithm  as  indicated 
in  the  above  examples  lead  us  to  believe  the  algorithm  will  become  a  popular 
method  for  solving  huge  positive  definite  Toeplitz  systems  of  equations. 
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